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Abstract 


A linearly implicit difference scheme for the space fractional coupled 
nonlinear Schrédinger equation is proposed. The resulting coefficient ma- 
trix of the discretized linear system consists of the sum of a complex scaled 
identity and a symmetric positive definite, diagonal-plus-Toeplitz, matrix. 
An efficient block Gauss-Seidel over-relaxation (BGSOR) method has been 
established to solve the discretized linear system. It is worth noting that the 
proposed method solves the linear equations without the need for any sys- 
tem solution, which is beneficial for reducing computational cost and mem- 
ory requirements. Theoretical analysis implies that the BGSOR method is 
convergent under a suitable condition. Moreover, an appropriate approach 
to compute the optimal parameter in the BGSOR method is exploited. Fi- 
nally, the theoretical analysis is validated by some numerical experiments. 
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1 Introduction 


The Schrédinger equation is a crucial equation in quantum mechanics, a 
science that studies submicroscopic phenomena. It can arise from the Brow- 
nian path integral. In [6], the path integral method to the Lévy-a@ process 
was generalized, and the space fractional equations were derived. 

Consider the space fractional coupled nonlinear Schrédinger (CNLS) 
equations 


a,<a4<a2, O0<t<T. 


(1) 


un (lu? + dlv/2) w=, 
v+n({v|? + 4|ul?) v =0, 


Given the conditions of the initial boundary value as follows: 


u(x, 0) =uo(#), v(x, 0) = vo(2), ays 
0 


< 
v(ai,t) =u(a2,t) =0, v(ai,t) = v(a2,t) = 0, t 


a2 
< T, 

where 2 is the imaginary unit, € > 0, 7 > 0, 8 => O are some constants, 
and 1 <a < 2. In [5], the fractional Laplacian was designated as 


(—A)? u(a,t) = #7 (|9|* (ula, t))), 


in which # stands for the Fourier transform applied to the spatial variable x. 
Assuming that _.. D?u(a, t) and ,D¢,,u(, t) are the left and right Riemann— 
Liouville fractional derivatives of order a € R* given by 


1 o” 7 n—-l-a 
-oDzu(z,t) = race on | 2-9 u(y, t)dy, 
2D .u(z,t) = raya | (u— a)" “u(u, t)dp, 


respectively, the Riesz fractional derivative can be calculated as 


lea 


aime t) = —(—A)? u(a,t) = — [-coD2u(z, t) + 2D?.u(x, t)] ; 


igo 
2 cos 73° 


In general, analyzing and understanding the behavior of the exact solu- 
tions of the space fractional CNLS equations is so challenging. In recent years, 
some numerical methods have been proposed to solve the CNLS equations. 
The difference method [12, 13, 11], the Crank—Nickelson scheme [1], and the 
collocation method [2] have been presented to solve the CNLS equations. 

The discretization of the CNLS equations leads to the solution of the com- 
plex symmetric linear systems. The coefficient matrix consists of the sum of 
the symmetric positive definite, diagonal-plus-Toeplitz, matrix and the com- 
plex identity scaled matrix. Recently, Dai and Wu [4] developed a suited 2x2 
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linear system and employed the block Gauss-Seidel (BGS) iteration scheme 
to solve the resulting linear systems. Then they analyzed the convergence 
of the BGS scheme for the corresponding 2x2 linear system. In this work, 
we establish a fast block Gauss-Seidel over-relaxation (BGSOR) scheme for 
solving the two-by-two linear system that arises from the discretization of 
CNLS equations. Notably, the new method allows the corresponding sys- 
tems to be solved without the need to compute the inverse of the coefficient 
matrices. Moreover, it should be pointed out that the BGS method can be 
regarded as a special case of the new method when the relaxation parameter 
is set to be one. 

The arrangement of this work is as follows. In Section 2, the model 
problem will be studied, and a linearly implicit difference technique will be 
presented. Application, convergence theory, and finding the optimal param- 
eter for the BGSOR method are proposed in Section 3. Section 4 is devoted 
to giving some numerical examinations. In Section 5, we finally made some 
conclusions. 


2 Model problem and a linearly implicit difference 
scheme 


The domain 2 = (a1, a2) x (0, T) is divided into a uniform grid of mesh points 
(x,;,t,), where 


xj =a, + gh, h= 


and 


th, =ktT, T= 


At grid points, the values of functions u(z,t), v(a,t) are, respectively, denoted 
by uy =u(x;, tz), ve =v(x;,th), and Ue, 5 are the approximate solutions 
of (1). 

The following equation gives a discrete approximation to ae u(a,t) [10]: 


oF e(e (a 
ane" u(xj,th) = er aly a ha( (aj— ate)+>_ wr alert) +0(h7), 
1=0 
(2) 
where VU, = Teo) and {w¢} is defined as follows: 
2 
~(a a a ~(a Qa a 
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Ortigueira [7] proposed the following fractional central difference operator: 
co 
pu(e) = S> 9 u(w—Ih), 
l=—oo 
where 


gl = (-1)'T(a+) 
: (¢—14+ D($ +141) 


As stated in [7], the coefficient {go} satisfies 


x 
2sin | — 
an(2) 


When a > —1, the recursive relations for { gh} are as follows: 


(a) _ T(a+ 1) (a) _ _ a+l1 .(@) as 
Io —_ T2(a/2+1)’ g — {—15 ey 


2 fore) 


— S- ger, ceER. 


l=—oo 


Lemma 1. [10] Assume that u(x) € C°(R) and that its all derivatives of 
order up to 5 belong to L1(R). Then, it holds 


Afu(z) — O°u(x) 


+ O(h?). 
ia Aiz|2 O(h*) (3) 
From Lemma 1, it follows that 
a Afu av 1 = a(@ 
(A) Fu(ey, ty) = SEM 4 O(n) = Kr al(e.t4) + OF) 


Now, we consider the following numerical scheme for solving (1) [12]: 


Up US m . Ai Uk 1 
AE Eyal (SE) + olla + arr) 


2T + he = 2 
x apr a wn 1 = 
9 a ’ 
(4) 
yktl _ yk-l 7 m «eg WE + ae 1 
PE Ea (AP) + otto + rege) 
I=1 
okt + yk-l 
+ ; =0, 
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where 1 < 7 < m, 1 < k <mn-—1. Another scheme should be provided for 
the numerical solution at k = 1. We consider the following scheme for this 
purpose (see [3]): 


U— Us ue 
i ¥ g6® (1) 2 0)2\ gl _ 
—— +e 2G US” + (IP? + BYP?) 
yili_wyo m 
Y (a 
r2£—4 4 HSH) + oll? + BlUPP)VH =0, 


Ui — YP ih Ui + UP 
j jd oy (a) 
ee +2 yale (= 3 ) 


3 1 3 1 WU) + Up 
o( Slap? - Brae? + a( S19fP1? - civ?) J =o, 


J J 2 2 
yi rr V1 4 0 
j j a ae a(a) { 71 l 
a +e LAH 71) 
3 1 3 e+ 
; yh |? 4 ys wl — 
o( Zing ~ Sipe + a( Freep? - Zia?) ) 25 


The structure of the first and second difference equations in (4) is the same. 
Set 


k+1 _ k+1 k+1)T k+1 _ saa k+1)T 
U =|[% puaes Dn ] ’ b = [b; iat leg ] ? 


ET ; : ; ; 
=jar Ge =mr(IZP +BY), DE =diag(dy*,....dm""). 


So, at each time step, we need to solve the following systems of linear equa- 


tions: 
Abtlg etl @ pet l<k<n-l, (5) 
Betiyktl — oktl 1l<k<n-l, 


b 


where A*t! = T + D¥+! + oF and b**! is as follows: 


WE — SAI — aa 


= 
Us = Oe ile _ dktlayk-l 
1=1 
peti = 


ana. 0 ga 
Uk 1 ~ nda UF 1 = diigk-1 


Moreover, T is the Toeplitz matrix, which has the following structure: 
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gl)... ghd. gh? 


G0? BS Gn 

T=pu : (6) 
Gene Geng + Gy OY 
(a) a(a) (a) 


i) al)... gl) 
Im—1 Im—2 gy Io 


Also, it should be noted that B*+! and c*+! can be obtained by changing 
the roles of Y and Y in A**! and b**1, 


3 The BGSOR iteration method 


To establish the BGSOR iteration method, we need to give some preliminar- 
ies. Let us first consider the iterative solution of the linear equation 


AY =b, (7) 
in which A € C’* is a nonsingular complex symmetric matrix as follows: 
Ape Day 


where T is the symmetric positive definite (SPD) and Toeplitz matrix des- 
ignated in (6), D = diag(d,,do,...,de) with d; > 0,7 = 1,2,...,2@, is the 
diagonal matrix, U,b € C’. Let U = x+y and b = f + 19 be complex 
vectors, where y,z,p,q € R‘. So, the system can be rewritten as a particular 


form, namely, 
y = es (Y)-(F)=2 (8) 
~LQT)] \« ee 


where Q = D+. We are now in a position to design a new method for 
solving (8). 


To introduce the BGSOR iteration method, we consider the following 
decomposition for the coefficient matrix (8): 


A = (wR - @)— (EF —(1-w)D) =: H-YN, (9) 


-10 0 0 
e-(G7) #= (90): 


and w is a positive parameter, which is known as the relaxation parameter. 
Using the decomposition (9), the BGSOR iteration method is stated as 


where 


Mit) = NZ") +B, k=0,1,2,..., 
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where .@ and .¥ are defined as (9), and z“*) = (y“); x"). Note that y“) and 
x‘*) are two M-vectors that stand for the iterations. Also, z© is an arbitrary 
initial guess. Thereupon, the iterations take the following procedure: 


ae = : ((w— Dy +Qx® — Ff), (10) 


((w 1)x™) +g Oy): 


As can be seen, there is not any system solution in each iteration, and only 
two matrix-vector multiplication are needed. This can be very important 
because the new scheme requires insignificant computational efforts and just 
contains the matrix-vector multiplications. Furthermore, if w = 1, then the 
iteration scheme (10) reduces to 


yt) — Ox) — Ff, as 
x(t) — g_ Gy (e+), 


which is presented in [4] and known as the BGS iteration method. Therefore, 
the BGS iteration method is a special case of the BGSOR iteration method. 

Next, we investigate the convergence of the BGSOR method for solving 
(8), and then we obtain the optimal value of the relaxation parameter w. In 
the following, we recall a result that will be useful later. 


Lemma 2. [14] Suppose that the quadratic equation x? — pr + q = 0, where 
p and q are real numbers. Both roots of the equation are less than one in 
modulus if and only if |g] < 1 and |p| < 14+ q. 


Theorem 1. Consider A = D+ T +21 € R®™** as a matrix, where D and 
T are diagonal and Toeplitz SPD matrices, respectively. The necessary and 
sufficient condition for convergence of the BGSOR iteration method to the 
solution of (8) for any initial guess, is 

i a Hmax(Q) 


an ane 


where Umax(Q) is the largest eigenvalue of Q. 


Proof. Let \ be an eigenvalue of the iteration matrix Y = .@~1.V, and let 
x = [u;v] be the corresponding eigenvector. Without loss of generality, let 
A #0. Then 

(D—w&)(6T — (1— aQ))x = dx, 


equivalently, 


(1—w)u—- Qv =— Awu, (12) 
(w—1)v =A\(Qu+ wv). (13) 


We can derive from (12) and the positive definiteness of Q that 
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v=((A-—1lw+1)Q*u. (14) 
Substituting (14) into (13), gives 
—rQ?u = ((A—1)w + 1)?u. (15) 
This shows that if yw is an eigenvalue of Q, then 


Ne Ss (OS Geet? (16) 
== (Aw? + Qu(1 — w)A + (w — 1)?). (17) 


From (17), we get 


? (HS) pCa 6, (18) 


Ww? Ww 


Now it follows from Lemma 2 that |A| < 1 if and only if 
lw — 1] <u, 
[Qu — Quy — p?|7 < Qu? — Qw +1. 


It is straightforward to see that |w — 1| <w is equivalent to w > $. By some 
easy manipulations, we can observe, whenever 


(2u = 1)? > 12, (19) 
the second inequality holds. The inequality (19) is ensured, if 
\2u —1] > or |2w—1| < —p, 


equivalently, 
1- 1 
w< SF or w> te (20) 
Evidently, the first inequality of (20) cannot be true. On the other hand, 
holding the second inequality of (20) ensures w > 5, and then it completes 


the proof. 


In the following, we would like to find the optimal value of the relaxation 
parameter w, denoted by w*. To do so, w* should be computed to minimize 
the spectral radius of the iteration matrix of the BGSOR method, that is, 


p(G%.*)=arg min p(G,). 


is a tHmax(@} 


To compute the optimal value of w, we state and prove the next theorem. 


Theorem 2. Assume that the hypothesis of Theorem 1 are met. Then the 
optimal value of the relaxation parameter and the corresponding optimal 
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convergence factor in the BGSOR iteration method are as follows: 


w= 5 (1+ VIF FO), (21) 


_ Le. p(Q) ° 
Bigs aaa Ao) : 


and 


Proof. If A is an eigenvalue of the iteration matrix Y,, then A < 0 or A € 
C\R, according to (16). First, we consider the case \ < 0. So, there exists 
an eigenvalue y of Q such that (18) holds true. The discriminant of this 
quadratic equation is 


a= (awa) 1(2) 
Ww wW 


and the roots of (18) are as follows: 


Dj? = Daya? af 
A1,2(w) = Dye =e 5) 7 


From (16), we get 


(A—1)w+1=+pvV-d. (22) 
Set 


f(A) =Q — 1)w+1=0dA4+1-y, 
g(A) =a pwr. 


Clearly, the function f,, passes through the point (1,1), that is, f,(1) = 1 
and the slope of f,,(A) is w. Figure 1 displays the points of intersections of 
the functions f,,(A) and g(A) for an arbitrary value of w. This figure shows 
that by increasing w, the maximum of absolute values of the abscissas of the 
points of intersection of the functions f,,(A) and g(A), that is, max{A1, Az}, 
decrease, while f,,(A) gets tangent to g(A). In the tangent case, we have 
Ay = re, and it indicates that A = 0. From A = 0, it is straightforward to 
verify that 2 = 0 or 4w? —4w— pu? = 0. The case ps = 0 is impossible, because 
of the positive definiteness of Q. Thus, 4w? — 4w — yw? = 0. This quadratic 
equation has two roots, as follows: 


ws = 5 (1+ VIF). 


Due to the condition w > “4 , w— is not acceptable. So, we consider 


1 
wy = 5 (1+ VIFe), 
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and in this case, we have 


1 
Ay =A. =AG=zE— - 1. 
WW. 
Now suppose that w > w+. In this case, the roots of the quadratic equation 
(18) are complex and conjugate, which are as follows: 


—2u? +2w+p? | vA’ 


Ee we a e 
where 2 
a'=4(—) (ee) 
Ww 7) 
Then ‘ 
|Ai2| =1-——. 
7) 
By recalling that w > w4 and having in mind that w; > 1, we have 
1 1 
1—— <1--, 
Wi WwW 


and this shows that w+ is the best choice for w. On the other hand, the curve 
g(A) = £p(Q)V—A serves an upper bound for each curve as +p./—A, where 
0< yu < p(Q). Summarizing the above results, we see that 


2 
p(G%*) =min max |1 ries4 : ={ p(Q) ) 


Dn yp beige wr \14 /1+ PQ) 


where w* was considered as in (21). 


Remark 1. In Theorem 2, for computing w*, we need to compute p(q). One 
may use a few iterations of the power method to compute Amax(Q). On the 
other hand, because of positive definiteness of Q, we have 


p(Q) = Anau) = I|Qll2- 


So, we can compute ||Q||2 instead of p(Q). In practice, the normest command 
of MATLAB can be used to compute an estimation of ||Q]]2. 


4 Numerical experiments 


This section is devoted to numerical experiments to evaluate the effectiveness 
of the BGSOR iteration scheme for solving linear systems (8). The numerical 
results of the proposed method are compared with those of the GMRES [8, 9] 
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Figure 1: The graph of the functions f,,(A) and g()). 


and the BGS methods. In all the test problems, we use the restart version 
of GMRES with a restarting number 10. The initial guess is assumed to be 
a random vector, and iterations are terminated when 


Res = (alle - 10-9, 


~ [rolls 


where r, = Y — &z*) is the residual at the kth iteration or if the maximum 
number of iterations mazxit = 1000 is exceeded. The terms “IT” and “CPU” 
in the tables refer to the total number of iterations and the elapsed CPU time 
in seconds for convergence, respectively. We comment that five runs were 
performed for each test, and then the average of CPU times and iterations 
are reported (The average of the iteration numbers were rounded). For the 
BGSOR method, the optimal parameter is computed according to the rule 
(21). The numerical results were carried out under MATLAB-R2017 on a 
laptop running Windows 10 and an Intel (R) Core(TM) i5-8265U CPU @ 
1.60 GHz 8 GB. 


Example 1. Let 6 = 0. The system (1) is then decoupled and becomes 
wy + (—A)2u+ 2|ul?u = 0, 

when the initial value 
u(x, 0) = sech(a) - exp(2z), 


is applied. In this example, the original problem was truncated in [—20, 20]. 
Set u(—20,t) = u(20,t) = 0. For this problem, we choose the parameters 
€=1.3 and 7 = 1.2. 
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Table 1: The optimal parameters w* for BGSOR method with a = 1.3 and n = 4m at 
t = 2 for Example 1. 


ro 


800 1600 3200 6400 
w* 1.002 1.004 1.006 1.009 


We set m = 800, 1600, 3200,6400 and examine two values of a, a = 
1.3,1.6. When a = 1.3, we set n = 4m; otherwise, we choose n = 6m. The 
optimal values of the relaxation parameter in the BGSOR method for a = 1.3 
are given in Table 1, and the ones for a = 1.6 are given in Table 3. 

In Tables 2 and 4, we have listed the numerical results at t = 2. From 
these tables, we observe that the BGSOR method is superior to the examined 
methods in terms of both the iterations and the elapsed CPU times. 


Table 2: Numerical results with a = 1.3 and n = 4m at t = 2 for Example 1. 


Method 7 300 1600 3200 6400 
Boer cnt ee On naa oa 
Bee ane ae 6p ee ie 
GMRESUY) ceo G8 ae ae re 


Example 2. For the following coupled system with 0 F 0: 


4 (_A)S a Di bee| Danae 
ae ORS lal eel) G, AVS EO ESD, 
wy, +(—A)?v42 (\v|? + \v|?) v=0, 
(23) 
We will use 


u(x, 0) = sech(a + Do) -exp(tvoxr), v(a,0) = sech(x — Do) - exp(—2voz), 
u(—20,0) = u(20,0) = 0, v(—20,0) = v(20,0) = 0, 

(24) 
as the initial conditions. In this case, we choose the parameters Do = 1, 
vo = 2, € = 1.4, and 7 = 1.2. 

The discretization of the coupled system of (23) leads to the solution 
of the linear systems of equations of the form (5). We assume that these 
coefficient matrices are A and B. These matrices have the same structure. 
Tables 5 and 7 show the optimal values of the relaxation parameter of A and 
B in the BGSOR method for different values of a and m. 
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Table 3: The optimal parameters w* for the BGSOR method with a = 1.6 and n = 6m 
at t = 2 for Example 1. 


£ 800 1600 3200 6400 
w* 1.010 1.022 1.050 1.108 
Table 4: Numerical results with a = 1.6 and n = 6m at t = 2 for Example 1. 


Method L 800 1600 3200 6400 
Boeor ane ae nae ali ane 
BoP opt cop nie oof ve 
GME UO) ce ve nas NOE soni 


In Tables 6 and 8, we report the results for the BGSOR, BGS, and GM- 
RES(10) iterative methods at t = 2. These results clearly show that the 
BGSOR method leads to a faster overall convergence time than the other ex- 
amined methods. Besides, the BGSOR method gets less iteration numbers. 


Table 5: The optimal parameters w* of A and B for the BGSOR method with a = 1.3 
and n = 4m at t = 2 for Example 2. 


U 800 1600 3200 6400 
w*(A) | 1.002 1.004 1.006 1.008 
w*(B) | 1.002 1.004 1.006 1.008 


Method L 800 1600 3200 6400 

A B A B A B A B 

IT 5 5 5 5 5 5 5 5 
ahaa: CPU 0.013 0.010 0.052 0.023 0.173 0.145 0.938 0.841 

BGs IT 5 5 6 6 6 6 7 7 
CPU 0.020 0.014 0.069 0.064 0.213 0.228 1.641 1.145 

; IT 6 6 7 7 7 7 8 8 
CMRES CPU 0.064 0.017 0.093 0.049 0.155 0.139 2.812 1.377 


5 Conclusion 


In this paper, the BGSOR scheme has been presented to solve the com- 
plex symmetric linear systems deriving from the discretization of the space 
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Table 7: The optimal parameters w* of A and B for the BGSOR method with a = 1.6 
and n = 6m at t = 2 for Example 2. 


L 800 1600 3200 6400 
w*(A) | 1.010 1.022 1.050 1.122 
w*(B) | 1.010 1.022 1.050 1.122 


Method L 800 1600 3200 6400 

A B A B A B A B 

IT 6 6 7 7 9 9 10 10 
Boece CPU 0.021 0.017 0.071 0.069 0.346 0.248 1.941 2.003 

BGS IT i ‘4 10 10 15 15 35 35 
CPU 0.025 0.020 0.106 0.112 0.607 0.592 6.832 5.483 

: IT 8 8 9 9 11 11 13 13 
ee) CPU 0.088 0.061 0.093 0.082 0.448 (0.412 3.376 3.251 


fractional CNLS equation. We have analyzed the convergence theory of the 
method, and we have shown that the method is convergent under a suitable 
condition. The optimal value of the relaxation parameter and the rate of 
convergence factor for the BGSOR method were also provided. Our results 
have verified that the BGSOR method performs better than some existing 
methods. 
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